Using state machines and callbacks to model friction and limit stops
Modeling friction is always an annoying task, requiring extra characterization of your system and convincing modeling tools to operate through all the nonlinearities of the model.
Here we'll make a state machine to think about the different states of the model and then use the VectorContinuousCallback feature of the Julia DifferentialEquations.jl package.
Consider a sliding mass system in which an external force is applied to the mass and the mass is subject to static and kinetic friction, and bounded by limit stops which abruptly stop its motion.
We can identify five different states for the system:
LEFT_STOP: mass against left stop, motion resisted by static friction and stopSLIDING_RIGHT: mass moving, being resisted by dynamic frictionSTOPPED: mass not against either stop, motion resisted by static frictionSLIDING_LEFT: likeSLIDING_RIGHTbut in the other directionRIGHT_STOP: likeLEFT_STOPbut at the other end.
There are ten transitions between states, organized in eight lines:
LEFT_STOPtoSLIDING_RIGHTif force (on the mass) is greater than static frictionSLIDING_RIGHTtoSTOPPEDif speed drops to zeroSLIDING_RIGHTtoRIGHT_STOPorSLIDING_LEFT(depending on COR) if position reaches right stopSTOPPEDtoSLIDING_RIGHTif force is greater than static frictionSTOPPEDtoSLIDING_LEFTif force is greater than static friction (but in the other direction)SLIDING_LEFTtoSTOPPEDif speed drops to zeroSLIDING_LEFTtoLEFT_STOPorSLIDING_RIGHTif position reaches left stopRIGHT_STOPtoSLIDING_LEFTif force is greater than static friction (but to the left)
Let's enumerate our states and make a default parameter object using some tools from Parameters.jl.
xxxxxxxxxx State LEFT_STOP SLIDING_RIGHT STOPPED SLIDING_LEFT RIGHT_STOPParamsxxxxxxxxxx mutable struct Params M::Float64 = 1.0 # mass of sliding object in kg ga::Float64 = 9.81 # gravitational constant, m/s^2 μs::Float64 = 0.8 # coefficient of static friction μk::Float64 = 0.5 # coefficient of dynammic friction cor::Float64 = 0.3 # coefficient of restitution, 0=inelastic, 1=elastic xls::Float64 = -0.5 # location of left stop in m xrs::Float64 = 0.5 # location of right stop in m state::State = STOPPED # state of systemendLet's write a state-aware differential equation, and write a separate function for calculating the net force which can be used by both the ODE and by the callbacks for interpolating to the state change points.
sbode (generic function with 1 method)xxxxxxxxxxfunction sbode(u,p,t) v,x = u # unpacking state variables p # unpack parameters by name into function scope #netforce = calcnetforce(p, force_ext, t) if state==LEFT_STOP dv = 0.0 # not moving dx = 0.0 # not moving elseif state==SLIDING_RIGHT dv = 1/M*(force_ext(t)-M*ga*μk) dx = v elseif state==STOPPED dv = 0.0 dx = 0.0 elseif state==SLIDING_LEFT dv = 1/M*(force_ext(t)+M*ga*μk) dx = v else # state==RIGHT_STOP dv = 0.0 dx = 0.0 end [dv, dx]end The tests we use to move from one state to the other are defined in FricSMCondx, while the corresponding actions are defined in FricSMaffect!.
FricSMCondx (generic function with 1 method)x
function FricSMCondx(out, u, t, integ) integ.p v,x = u out[1] = (state==LEFT_STOP)*(force_ext(t)-M*ga*μs) out[2] = (state==SLIDING_RIGHT)*(-v) out[3] = (state==SLIDING_RIGHT)*(x-xrs) out[4] = (state==STOPPED)*(force_ext(t)-M*ga*μs) out[5] = (state==STOPPED)*(-force_ext(t)-M*ga*μs) out[6] = (state==SLIDING_LEFT)*v out[7] = (state==SLIDING_LEFT)*(-(x-xls)) out[8] = (state==RIGHT_STOP)*(-(force_ext(t)+M*ga*μs))endFricSMaffect! (generic function with 1 method)x
function FricSMaffect!(integ, idx) integ.p print("affect,",integ.t,",",idx) if idx==1 integ.p.state = SLIDING_RIGHT elseif idx==2 integ.p.state = STOPPED integ.u[1] = 0.0 elseif idx==3 if cor==0.0 integ.p.state = RIGHT_STOP integ.u[1] = 0.0 else integ.p.state = SLIDING_LEFT integ.u[1] = -cor*integ.u[1] end elseif idx==4 integ.p.state = SLIDING_RIGHT elseif idx==5 integ.p.state = SLIDING_LEFT elseif idx==6 integ.p.state = STOPPED integ.u[1] = 0.0 elseif idx==7 if cor==0.0 integ.p.state = LEFT_STOP integ.u[1] = 0.0 else integ.p.state = SLIDING_RIGHT integ.u[1] = -cor*integ.u[1] end elseif idx==8 integ.p.state = SLIDING_LEFT endendx
cb = VectorContinuousCallback(FricSMCondx, FricSMaffect!, nothing, 8); Create a forcing function.
force_ext (generic function with 1 method)xxxxxxxxxxforce_ext(t) = t<25 ? t/2*sin(2π*t/2) : (25-t/2)*sin(2π*t/2)Set up a sim function which sets initial conditions, timespan, parameters, and simulation options. We set dtmax = 0.01 because without it the solver skips forward too quickly.
sim (generic function with 1 method)x
function sim() u0 = [0.0, 0.0] tspan = (0.0, 50.0) p = Params() prob = ODEProblem(sbode, u0, tspan, p) sol = solve(prob, callback = cb, dtmax=0.01)endxxxxxxxxxxsol = sim();Plot the results! Here we can see the system being pushed back and forth by an amplitude-modulated sinusoidal force.